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We present a fast one dimensional hybrid method to efficiently simulate extensive air showers up 
to the highest observed energies. Based on precalculated pion showers and a bootstrap technique, 
our method predicts the average shower profile, the number of muons at detector level above several 
energy thresholds as well as the fluctuations of the electromagnetic and hadronic components of the 
shower. We study the main characteristics of proton-induced air showers up to ultra-high energy, 
comparing the predictions of three different hadronic interaction models: S1BYLL 1.7, S1BYLL 
2.1 and QGSjet98. The influence of the hadronic interaction models on the shower evolution, in 
particular the elongation rate, is discussed and the applicability of analytical approximations is 
investigated. 



PACS numbers: 13.85.-t,96.40.Pq,96.40.-z 



I. INTRODUCTION 



Extensive air showers (EAS) generated by cosmic rays 
in the Earth's atmosphere are the only way to study cos- 
mic rays of energies above 10 15 eV. At lower energies the 
cosmic ray spectrum and composition are studied in ex- 
periments that measure directly the charge and energy 
of the primary particle. The analysis of air shower data 
relies on simulations that use the current knowledge of 
hadronic interactions to predict the observable shower 
parameters. With increasing cosmic ray energy, this task 
becomes more difficult as the gap between the shower en- 
ergy and the energy range studied in accelerator exper- 
iments increases and the hadronic interaction properties 
have to be extrapolated over a wide range. The difficul- 
ties are also related to the fact that particles produced in 
the forward region of the interaction are not registered in 
collider experiments, while they are responsible for most 
of the shower characteristics. Last but not least, the at- 
mospheric targets are light nuclei which have not been 
studied in collider experiments. 

Air shower experiments are either ground arrays of de- 
tectors that trigger in coincidence when the shower passes 
through them, or optical detectors that observe the lon- 
gitudinal development of EAS. Both types of instruments 
are sometimes supplemented by shielded or underground 
detectors that observe the muon component of the show- 
ers. The most commonly observed EAS parameters are 
the number of charged particles at ground level for the 
shower arrays, or at shower maximum (S'max) for the op- 
tical detectors; the depth of shower maximum (X max ) 
itself, and the number of muons (N^) above different en- 
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ergy thresholds. The combination of these and occasion- 
ally additional shower features, calculated in simulations 
with a particular hadronic model, is used as the basis for 
the determination of the energy and mass of the primary 
particle. Reviews of air shower experiments and observed 
features are given, for example, in ||]. 

At the end of the cosmic ray spectrum, at energies 
above 10 19 eV, air shower simulation becomes a very dif- 
ficult problem technically. The number of charged parti- 
cles that have to be followed in the Monte Carlo scheme 
is proportional to the shower energy. For example, high- 
est energy cosmic ray showers |3[ Q can have more than 
10 11 charged particles at X max . As a consequence the 
direct simulation of the shower following each individual 
particle becomes practically impossible, especially when 
a large number of showers has to be simulated. 

The widely used solution to the problem of having to 
deal with an excessively large number of shower particles 
is the simulation of EAS using the thinning technique f| . 
This method is extremely useful to estimate detectable 
signals and to compute average values of the observables 
1^, 0. The thinning procedure follows only a subset of 
the shower particles below a certain energy threshold, 
assigning weights to them so that the average number 
of particles at the ground is correctly reproduced. Due 
to this, artificial fluctuations are introduced even when 
small energy thresholds are used. Various methods of re- 
ducing artificial fluctuations have been proposed recently 
(e.g. j^, ^|) optimizing the compromise between time- 
consuming simulations and fluctuation-enhancing thin- 
ning. 

In this work we present a hybrid method of simulating 
the longitudinal profile of extensive air showers. It is a 
fast, one dimensional calculation which provides predic- 
tions for the total number of charged particles and muons 
along the shower axis. The method allows the collection 
of sufficiently high Monte Carlo statistics without losing 
information about shower fluctuations. 

In general, hybrid calculations are based on the idea 
to follow the development of air showers in detail above 
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a certain energy threshold and to replace subthreshold 
particles by a simplified and efficient approximation of 
the subshowers initiated by them. Many hybrid calcu- 
lations use the Monte Carlo-generated high-energy sec- 
ondary particles of the first few interactions of a cosmic 
ray in the atmosphere as initial distribution, and then 
calculate the particle densities observed at detector level 
by solving the corresponding transport equations (see, 
for example, @ 0, g g @). 

Here we follow the approach of Gaisser et al. |15[ and 
treat the subthreshold particles with a library of shower 
profiles based on presimulated pion-initiated showers. 
This idea can be combined in a bootstrap procedure 
to extend the shower library to high energy. The novelty 
of this work is that we extend the method of [ll| by 
accounting for fluctuations in the subshowers generated 
with the shower library, and also calculate the number of 
muons at detector level above several energy thresholds. 

Showers simulated in this way can be used as input 
to simulations for experiments measuring the longitu- 
dinal development of the shower such as HiRes |L7) , 
the fluorescence detector of the Pierre Auger Observa- 
tory and future experiments such as EUSO |L9] , 
OWL/AirWatch |20) and the Telescope Array Be- 
sides this, as will become clear later, hybrid simulations 
are very helpful for comparing shower parameters pre- 
dicted by different hadronic interaction models and to 
aid the interpretation of the experimental results in this 
way. 

This paper is structured as follows: In Section II we 
describe the hybrid method and the parametrizations 
of the presimulated showers. We demonstrate the self- 
consistency of the method by comparing showers simu- 
lated directly with predictions from the hybrid calcula- 
tion. In Section III we apply the hybrid method to proton 
induced showers at fixed energy. We give the average val- 
ues and distributions of X max , S m&x , and obtained for 
different hadronic models and discuss how the differences 
are related to the simulation of high-energy multiparti- 
cle production. In addition the elongation rate theorem 
p2| , [2^ ] is discussed in terms of the different hadronic 
interaction models and their influence on the position of 
the shower maximum. Where available, we compare our 
predictions to calculations performed with the CORSIKA 
code Q which uses the thinning approach. Section IV 
summarizes our results and concludes the paper. 



II. THE HYBRID METHOD 

The hybrid method used in this work consists of cal- 
culating shower observables by a direct simulation of the 
initial part of the shower, tracking all particles of energy 
> fE, where E is the primary energy and / is an appro- 
priate fraction of it (in the following we use / = 0.01). 
Then presimulated showers for all subthreshold particles 
are superimposed after their first interaction point is sim- 
ulated. The subshowers are described with parametriza- 



tions that give the correct average behavior and at the 
same time describe the fluctuations in shower develop- 
ment. The method is extended recursively to higher en- 
ergies by bootstrapping. The results obtained at any 
primary energy can then be used for the simulation of 
showers at higher energy. 

It is well known that the fluctuations in shower prop- 
erties are dominated by fluctuations in the earliest 
and most energetic part of the cascade. We however 
parametrize both the average behavior and the shower 
fluctuations starting at 10 GeV. In this way we can use 
the hybrid method at relatively low energy of 100 - 1000 
TeV, where the results can be compared to those of direct 
(fully simulated) shower calculations. 

We build a library of presimulated showers by inject- 
ing pions of fixed energy E n , at fixed zenith angle 9 and 
depth X measured along the shower axis. The atmo- 
spheric density adopted here corresponds to Shibata's fit 
of the US Standard Atmosphere 1 25 , E6| , very similar to 



Linsley's parametrization. We limit the injection zenith 
angles to 9 < 45° since mainly showers in this angular 
range have been used for studies of the cosmic ray energy 
spectrum at the highest energies. 

Nucleon initiated showers are not presimulated. Nu- 
cleons are followed explicitly in the Monte Carlo down 
to the energy threshold for particle production. A sub- 
shower initiated by a kaon is assumed to be similar to 
one initiated by a pion of the same energy but with a dif- 
ferent first interaction point, which is sampled from the 
corresponding interaction length distribution. This ap- 
proximation is not expected to affect significantly our fi- 
nal results, the main reasons being the similarity between 
pion and kaon induced showers at high energy combined 
with the fact that the main contribution to shower de- 
velopment in this method comes from the highest energy 
particles that enter the parametrizations. Unstable par- 
ticles, including 7r°, i], A, £ and are allowed to interact 
or decay in the code. The interaction of these particles 
becomes important at the highest energies and account- 
ing for them can influence the average values of some 
observables. 

Photon and electron/positron induced cascades are 
treated with a full screening electromagnetic Monte Carlo 
in combination with a modified Greisen parametrization. 
The electromagnetic branch of the Monte Carlo includes 
photoproduction of hadrons. For energies above 1 EeV, 
the Landau-Pomeranchuk-Migdal (LPM) effect |27|, |f, 
p^ , p0{ is taken into account using an implementation by 
Vankov fl3~l| . The influence of the geomagnetic field on 
the cascade development Q is neglected. 

We have simulated primary pions of energies between 
10 GeV and 3 EeV with a step in energy of half a decade, 
interacting at fixed atmospheric depths Xq = 5, 50, 100, 
200, 500 and 800 g/cm 2 . For each pion energy, injection 
zenith angle and depth (i.e. a single entry in the library) 
we simulate 10,000 showers (5,000 at high energy) and 
record -X" max , .Smax, the longitudinal shower profile, and 
the number of muons above the threshold energies of 0.3, 
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FIG. 1: The correlation between lg X max and lg .S'max for 5,000 
pion induced showers at primary energy 3 x 10 18 eV initiated 
at Xq = 5 g/cm 2 and zenith angle 8 = 45° . 



1, 3, 10 and 30 GeV both at sea level and at a depth 
of 400 g/cm 2 above sea level measured along the shower 
axis. These values are used to produce distributions of 
showers in X max and S ma x, the correlations between them 
and distributions of the number of muons at sea level. 
The whole procedure of generating a library has to be 
carried out for each of t he int eraction models we adopt 
in this work (see section [II A). 



Fig. |j] shows an example of the correlation between 
X m &x and S'max- The plot contains 5,000 simulated pion 
showers of energy 3.16 x 10 18 eV initiated at atmospheric 
depth of Xo—5 g/cm 2 and zenith angle 9 — 45°. Correla- 
tions similar to these are produced for each entry in the 
library. Their correct representation is crucial for the 
successful modeling of shower fluctuations. 

Although it is unlikely to produce a high energy pion 
deep in the atmosphere, we also calculate their interac- 
tions at depths as large as 500 and 800 g/cm 2 to obtain an 
accurate description of the muon numbers at sea level and 
a better description of the late developing electromag- 
netic showers. For showers initiated after 500 g/cm 2 the 
atmosphere has been artificially extended beyond ground 
level. The distributions of muons are easily extended to 
other depths (corresponding to the observation level of 
different experiments) by extrapolation. For this task we 
use the slope of the muon longitudinal profile between 
sea level and a slant depth of 400 g/cm 2 above sea level. 

The longitudinal development of subthreshold meson 
induced showers is parametrized using a slightly mod- 
ified version of the well-known Gaisser-Hillas function 
that gives the number of charged particles at atmospheric 
depth X, pi: 



X-X 
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x exp 



(X - X n 



X(X) 



(1) 



Here X(X) = X + bX + cX 2 where Ao, b and c are treated 
as free parameters. Xq is the depth at which the first 
interaction occurs. The parameters b and c are assumed 
to be the same for all showers initiated at a given depth, 
angle, and energy. They are determined by fitting the 
mean shower profile of the parametrized showers to that 
obtained from the simulated shower profiles. 

The innovative approach of our method is that instead 
of using the average values of X m&K and S'max to gen- 
erate subthreshold meson showers of a certain energy, 
we sample their values (as well as the number of muons) 
from their corresponding presimulated distributions, tak- 
ing into account the correlation between them (Fig. [l]). 
This procedure accounts for the fluctuations in the sub- 
shower development. A technical remark is that we sam- 
ple the observables directly from their precalculated his- 
tograms, i.e. we do not assume any functional form for 
the distribution. In this way our code is very flexible - 
it allows the study of hadronic models that predict dis- 
tributions of observables not easily fitted by analytical 
functions. 

We sample meson subshowers at a zenith angle, depth 
and/or primary energy different from those we have pres- 
imulated by interpolating between the relevant parame- 
ters of the shower development (X max , S'max, Xq, 6, c, N^), 
corresponding to presimulated entries in the library 
which are adjacent in angle, energy, and depth to the 
subshower we want to describe. 
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FIG. 2: Energy distribution of pions in showers initiated by 
primary protons at E = 10 15 eV. The dashed curve is the 
energy distribution of the pions actually treated in our hy- 
brid simulation procedure using a hybrid energy threshold of 
10 13 eV. The solid curve shows the energy distribution of pi- 
ons which are explicitly tracked in a direct simulation. Pions 
which decay are not shown. 
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In Fig. g we plot the energy distribution of pions actu- 
ally treated in our hybrid simulation procedure (dashed 
line). For comparison we also show the energy distribu- 
tion of pions that must be explicitly tracked in a direct 
simulation (solid line) at the same energy. The compari- 
son is made for E — 10 15 eV proton induced showers and 
for the nominal energy threshold of 0.01 E. In the hybrid 
approach only the interactions of pions above 0.01 E are 
directly simulated and all lower energy pions are replaced 
with parametrizations. For a primary energy 10 15 eV we 
typically treat 1 out of 10 pions of energy ~ 3 x 10 10 eV as 
can be obtained from the figure. This explains the saving 
in CPU time achieved with the hybrid code with respect 
to the direct simulation - a factor about 7 for the par- 
ticular energy shown here. This factor rapidly increases 
with energy. Already at E = 10 16 eV the hybrid calcu- 
lation is about 25 times faster than a direct simulation 
(Table |). For applications which do not depend on the 
number of muons this factor increases even further. At 
10 15 eV about 25% of the CPU time is spent on correctly 
tracking the numerous muons that are decay products of 
charged pions and kaons which do not initiate hadronic 
showers, and hence do not enter the parametrizations. 
With increasing energy the number of mesons which de- 
cay at an energy above the hybrid threshold decreases, 
and the number of muons which have to be simulated 
explicitly becomes negligible. 

To ensure the consistency of our simulation approach, 
we have compared full simulations of pion showers to hy- 
brid simulations for the same initial energy and depth 
using several energy thresholds. We find a very good 
agreement between the average values of the different ob- 
servables and their fluctuations in the direct and hybrid 
simulations. Table | compares the direct simulations and 
the hybrid method for 5,000 vertical pion showers with 
fixed first interaction point at Ao=5 g/cm 2 , energy 10 16 
eV, and for the different hadronic models. It is very 
important to note that the differences between the two 
methods of calculation are much smaller than those intro- 
duced by the different hadronic interaction models, i.e. 
by using the hybrid approach we do not lose sensitivity 
to the models we are considering. 

In Fig. U we plot the distribution of the number of 
muons with energy above 0.3 GeV at sea level for vertical 
pion-induced showers of energy 10 16 eV. We compare the 
direct simulation to the results of the hybrid approxima- 
tion. Panel (a) shows this comparison for QGSjet98 and 
panels (b) and (c) are for SIBYLL 2.1 and SIBYLL 1.7, 
respectively. The relative differences in the average num- 
ber of muons are less than 0.5% for all hadronic interac- 
tion models. The same comparison for showers generated 
by primary pions with incident zenith angle of 45° shows 
larger differences between direct and hybrid simulations, 
but they are smaller than 2%. We believe these relatively 
small errors come mostly from the representation of the 
intrinsic fluctuations in the shower development and from 
the interpolation in energy and atmospheric depth that 
the code performs. Due to the bootstrap-like calculation 
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FIG. 3: Shower distribution in number of muons of energy 
above 0.3 GeV at sea level. Results are shown for 5,000 ver- 
tical pion-initiated showers of energy E = 10 16 eV, at fixed 
interaction point Xq = 5 g/cm 2 for different hadronic mod- 
els. The solid line represents fully simulated showers while 
the dotted line shows hybridly simulated showers with meson 
energy threshold _E/100. 



of the high-energy part of the library this error increases 
slowly with energy. Comparisons with direct simulations 
at 10 18 eV show that the deviations in the average num- 
ber of muons is typically 4%. Even for QGSjet98, which 
predicts the largest fluctuations at this energy, we obtain 
a very good description of the distributions. Their width 
is reproduced with an error smaller than 3%. The hybrid 
code always tends to underestimate the number of muons 
and its fluctuations. 

Our results also show a remarkable stability under 
changes of the energy threshold, from which we conclude 
that the primary to threshold energy ratio we have used 
(-E-thr = E/100) is sufficient to achieve a very good de- 
scription of the average values and fluctuations of ob- 
servables in nucleon and pion initiated showers. Using a 
threshold for mesons and for the electromagnetic compo- 
nent fixed to EX = = ^/10, we still obtain a good 
agreement for the average values. However we might not 
correctly include some of the extreme fluctuations that 
are possible in the early development of the showers. 
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TABLE I: Average values of different observables and standard deviation of their distributions obtained by direct and hybrid 
simulations of 5,000 vertical pion showers with fixed interaction point Xq=5 g/cm 2 , and primary energy E — 10 16 eV. The 
predictions of SIBYLL 1.7, SIBYLL 2.1 and QGSjet98 are presented. The energy threshold in the hybrid calculation is 
0.01 E = 10 14 eV. 



Model 


SIBYLL 1.7 


SIBYLL 2.1 


QGSjet98 






A-JLL fciU L lly UI 1U- 


Tiircif f T-TT7"nT"i/"T 
LJ1L cL L IlVUllU 


<X max ) [g/cm 2 ] 
o (A max ) [g/cm 2 ] 


603 602 
49 50 


587 586 
51 49 


574 576 
55 56 


(5 max )/S [GeV- 1 ] 
a (S m ^/E) [GeV" 1 ] 


0.75 0.76 
6.8 x 10~ 2 6.8 x 10~ 2 


0.75 0.75 
6.3 x 10~ 2 6.2 x 10~ 2 


0.75 0.75 
6.5 x 10~ 2 6.5 x 10 -2 


<iV M > (> 0.3 GeV) 


5.39 x 10 4 5.41 x 10 4 
1.79 x 10 4 1.81 x 10 4 


6.10 x 10 4 6.13 x 10 4 
1.86 x 10 4 1.87 x 10 4 


6.87 x 10 4 6.91 x 10 4 
2.25 x 10 4 2.28 x 10 4 


CPU Time [min]* 


935 33 


1091 41 


1398 79 



*A11 CPU times illustrated in this work refer to a 1 GHz AMD Athlon processor. 



III. APPLICATIONS 



Hadronic Interaction Models 



In this section we apply the hybrid approach described 
above to simulate proton-initiated showers at fixed en- 
ergy. These showers show most transparently the influ- 
ence of the hadronic interaction model on air shower ob- 
servables. In a forthcoming paper we will calculate pre- 
dictions for a realistic cosmic ray spectrum with a mixed 
cosmic ray composition consisting of protons and nuclei 



In the following we consider the hadronic interaction 
models SIBYLL, and QGSjet. We have created libraries 
for the model versions SIBYLL 1.7 @, SIBYLL 2.1 
{H |H and QGSjet98 [||. QGSjet and SIBYLL are suf- 
ficiently different to illustrate various important points 
of how properties of hadronic interactions are reflected 
in shower observables. In addition they are commonly 
used for the analysis of air shower measurements. In 
discussing the models we will focus on QGSjet98 and 
SIBYLL 2.1 and show SIBYLL 1.7 predictions only for 
reference purposes because many air shower data have 
been already analyzed with this model. 



QGSjet98 and SIBYLL 2.1 were shown to describe well 
collider data up to the highest energies available so far 
(see for instance J39|). However, already the extrapo- 
lation to the unmeasured parts of the phase space is 
different. These differences are amplified by going from 
proton-proton/antiproton to proton-air collisions. 
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SIBYLL 2.1 shows a considerable improvement with 
respect to version 1.7 in describing the measurements 
of hadronic interactions at collider energies. The im- 
portant changes in SIBYLL are the implementation of 
new parton densities and parton saturation, a new model 
for diffraction dissociation, and an energy-dependent soft 
component p6| . Nevertheless, at the highest cosmic ray 
energies, its predictions are similar to those of SIBYLL 
1.7. On the other hand, QGSjet98 predicts a high-energy 
extrapolation which is strikingly different from that of 
SIBYLL. 



FIG. 4: Inelastic proton-proton and pion-proton cross sec- 
tions as predicted by QGSjet98 and SIBYLL. 

One of the key features of the hadronic interaction 
models is their prediction on hadron-air cross sections. 
The proton-air cross section determines the height of the 
first interaction in the atmosphere. However, it should 
be emphasized that the pion- and kaon-air cross sec- 
tions are also very important for the shower development. 
Fig. [| shows the model cross sections for proton- and 
pion-proton collisions which are the input for the calcu- 
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lation of hadron-air cross sections. In both models free 
parameters are adjusted to fit the measured pp and pp 
cross sections which cover the energy range from the low 
end up to £q a b ~ 1.7 x 10 15 eV, i.e. Tevatron center- 
of-mass energy of y/s = 1800 GeV. The model predic- 
tions for the pion-proton cross section diverge already at 
much lower energy. The experimental restrictions here 
are much smaller since the pion-proton cross section is 
experimentally known only up to -Eiab = 4 x 10 11 eV. 
The difference in the high-energy extrapolation of the 
models arises from different assumptions on the spatial 
distribution of partons in protons and pions. Both mod- 
els implement the eikonal approximation but differ in 
many technical details such as the treatment of inelas- 
tic diffraction. In the following we discuss only the most 
basic version of the eikonal model as it is sufficient for 
explaining the important differences. 

In the eikonal model the inelastic cross section is given 
by 

cr ine l = y > d 2 6(l-exp{-2x s (s,fo)-2 X ? i (s,6)}) , (2) 

where the eikonal function is written as the sum of soft 
and hard contributions, Xs and \h- The two-dimensional 
impact parameter of the collision and the squared center- 
of-mass energy of the collision are denoted by b and s. 
At high energy one has Xh 3> Xs and the inelastic cross 
section is dominated by Xh, written as 

Xh (s,b) = ^* Q cn(p? toS ,s)A(s,b), J d 2 bA(s,b) = 1 . 

(3) 

The normalized profile function A(s, b) describes the dis- 
tribution of partons in the plane transverse to the col- 
lision axis. The minijet cross section ctqcd depends on 
the collision energy and the transverse momentum cutoff, 
pcutoff^ needed to restrict the calculation to the perturba- 
tive region. For a given energy dependence of the minijet 
cross section, only the profile function A(s, b) determines 
the inelastic cross section and its energy dependence. 

Qualitatively, QGSjet is a model which assumes a 
Gaussian profile function EOj 



A(s,b) 



1 



irR 2 



exp 




(4) 



with R being a parameter. The SIBYLL model is based 
on the Fourier transform of the electromagnetic form fac- 
tor, assuming that the distribution of gluons in a hadron 
is similar to that of the quarks. The corresponding pro- 
file function is energy-independent and is, for example, 
for proton-proton scattering |3q] 



A (b) = —{v\b\fKMb\) , 



(5) 



where K 3 denotes the modified Bessel function of the 
third kind and i/aO.7-1 GeV" 1 . 



For all \b\ < b s with Xh(s, b s ) 3> 1 the saturation limit 
is reached. From (Q) it follows that any further increase 
of the minijet cross section would not change the con- 
tribution to the cross section integral from the impact 
parameter region |6| < b s . This allows us to give a rough 
estimate of the energy dependence of the inelastic cross 
section at very high energy. For a QCD cross section 
dependence of ctqcd ~ s A , as is expected within pertur- 
bative QCD PHI , one gets for a Gaussian profile 



and at high energy 



CTincl 



6 2 ~i? 2 Alns (6) 



d 2 b 6(b s - b) = 7ri? 2 Alns (7) 



For R being energy-independent the cross section will rise 
only logarithmically with the collision energy. However, 
the parameter R itself depends on the collision energy 
through a convolution with the parton momentum frac- 
tions, R 2 w i?Q + 4o^ ff In s and a' cS w 0.11 GeV -2 . Hence 
the QGSjet cross section exhibits a faster than In s rise 



CTincl 



4-7T Attgg- In 2 s 



(8) 



The cross section limit for SIBYLL can be derived in the 
same way from Eq. (||) 



CTincl 



A 2 

7rc^- In . 



(9) 



where the coefficient c rs 2.5 was found numerically. 

Both cross sections satisfy the Froissart bound and ex- 
hibit a In 2 s energy dependence. However, the numer- 
ical factors are different. Assuming A w 0.25, then 
4ira' eS A w 0.13 mb and c7rA 2 /^ 2 w 0.2 mb, which ex- 
plains the faster increase of the inelastic cross section in 
the SIBYLL model. A larger power of A w 0.4, as im- 
plied by data from the HERA collider [[l2| , even amplifies 
the model differences. However, the difference between 
the model predictions is smaller than expected from the 
arguments above, the reason being a somewhat smaller 
minijet cross section assumed in SIBYLL as compared 
to QGSjet. The saturation of the parton densities im- 
plemented in SIBYLL tames their rapid growth at small 
parton momentum ]3(| . 

Information on the profile function can be derived by 
comparing the differential elastic cross sections measured 
at accelerators with model predictions |l3[ Q . The form 
factor approach describes current data reasonably well 
p5[ , whereas a Gaussian profile shows large, systematic 
deviations and predicts a wrong curvature. Although 
currently available data clearly favor profile functions de- 
rived from electromagnetic form factors, it is not clear 
whether this approximation is still good at ultra-high en- 
ergy. 

For hadron-air collisions (Fig. ||) the relative uncer- 
tainty in the extrapolated cross sections is considerably 
smaller than that of proton-proton and pion-proton cross 
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sections. The geometrically large size of the target nu- 
cleus (mainly nitrogen or oxygen) dominates the interac- 
tion cross section. At the highest energy considered here 
the relative difference is less than 15%. 
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FIG. 5: Proton- and pion-air production cross section. The 
production cross section is denned as the cross section for all 
collisions in which at least one new particle is produced. It 
can be written as cr pro d = o"tot — c e i — o" qc i where o"tot is the 
total cross section and <7 e i and er qe i are the elastic and quasi- 
elastic cross sections respectively. 

The evolution of air showers in the atmosphere de- 
pends directly on how much energy is transferred in each 
hadron interaction into the electromagnetic component of 
the shower. It is common to describe this energy trans- 
fer in terms of the elasticity of the interaction. Fig. ^ 
shows the mean elasticity of proton- and pion-air inter- 
actions as predicted by QGSjet98 and SIBYLL 2.1. We 
define the elasticity of an inelastic interaction (including 
diffraction dissociation) as K c \ = E\ cat ±/ E pro j where -Ei ea d 
is the energy of the most energetic hadron with a long 
lifetime (i.e. proton, neutron, A, and charged pions and 
kaons) and E pro j is the energy of the projectile particle. 
SIBYLL 2.1 consistently predicts more elastic collisions 
than QGSjet98 with a relative difference of up to 17%. 
Assuming similar other characteristics of hadronic inter- 
actions, a model with large elasticity predicts air showers 
which develop deeper in the atmosphere. 

Other important aspects relevant to air showers are 
the predicted multiplicity of secondaries and the energy 
fraction carried by neutral 7r°'s, which are closely related 
to the elasticity. Neutral pions decay immediately into 
two photons and feed the electromagnetic component of 
the shower. At the highest energies some neutral pions 
also interact hadronically because of the enormous time 
dilation. On the other hand, the charged particle multi- 
plicity is a measure of how fast the initial energy is dis- 
sipated into many hadronic low-energy subshowers. It is 
also a good indicator for the muon multiplicity since the 
decaying charged pions are the primary source of muons. 
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FIG. 6: Mean elasticity in proton-air collisions as predicted 
by QGSjet98 and SIBYLL 2.1 (see text). 

Fig. |?j shows the mean charged particle multiplic- 
ity in proton- and pion-air collisions as calculated with 
QGSjet98 and SIBYLL 2.1. QGSjet98 predicts a power 
law-like increase of the number of secondary particles up 
to the highest energy. In contrast, the SIBYLL multi- 
plicity exhibits a logarithmic growth similar to In 2 s at 
high energy. In the energy region from 10 13 to about 
10 16 eV both models predict the same multiplicity in 
p-air collisions. However the pion-air multiplicities are 
significantly different at all energies. In SIBYLL differ- 
ent parton densities are used for pions and protons. The 
currently implemented parametrizations from Gliick et 
al. |^6|, predict fewer partons at low x in pions as 
compared to protons. The predicted secondary particle 
multiplicity is strikingly different at the highest energies. 
QGSjet98 predicts more than twice as many secondaries 
as SIBYLL. The multiplicity of neutral pions is closely 
linked to that of charged particles and hence shows qual- 
itatively the same behavior. 

The differences in multiplicity can again be qualita- 
tively understood by considering Eq. (|2j). The minijet 
cross section predicted by perturbative QCD describes 
the inclusive cross section of minijet pairs. It does 
not specify how many minijets are produced per single 
hadron-hadron collision. The mean minijet multiplicity, 
(rijct), can only be calculated after knowing the inelastic 
cross section 

(rijet) = 0-QCDMnel ■ (10) 

The larger multiplicity predicted by QGSjet stems both 
from the steeper energy dependence of its minijet cross 
section and from the more moderate energy dependence 
of its inelastic cross section. A detailed discussion of the 
relation between the minijet cross section and secondary 
particle multiplicity is given in p7| . 

Another difference is emphasized in the inset in Fig. [7]. 
At low energy (i.e. 100 to 1000 GeV lab. energy) the mul- 
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FIG. 7: Mean multiplicity of charged secondary particles 
produced in inelastic proton- and pion-air collisions. 



tiplicity predicted for proton-air collisions is up to 25% 
lower in SIBYLL than in QGSjet. Whereas this differ- 
ence is unimportant for electromagnetic shower variables, 
it becomes observable in the number of low-energy muons 
produced in the decay of charged pions and kaons. 
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FIG. 8: Mean energy fraction carried by neutral pions, elec- 
trons and photons in inelastic proton- and pion-air collisions. 



B. Shower Size and Depth of Maximum 

I max and S'max are two typical shower parameters 
measured by fluorescence and Cherenkov light detectors 
in several experiments. Knowing the shower energy, the 
mean depth of shower maximum and its fluctuations can 
be used to infer the primary cosmic ray composition. 

Fig. H shows the average value of Y max as a func- 
tion of primary energy for proton showers injected at a 
zenith angle 9 = 45°. The lines were produced averaging 
I max over 5,000 showers. The predictions of SIBYLL 1.7, 
SIBYLL 2.1 and QGSjet98 are shown. The first impor- 
tant feature is that SIBYLL 2.1 predicts smaller (Y max ) 
values than SIBYLL 1.7 by about 22 g/cm 2 from 10 14 
to 3 x 10 20 eV. The predictions of SIBYLL 2.1 are closer 
to the values produced by QGSjet98. In fact, at ener- 
gies below about 3 x 10 17 eV the difference is smaller 
than 10 g/cm 2 and it increases with energy up to a max- 
imum of 27 g/cm 2 at 3 x 10 20 eV. QGSjet98 predicts 
values of (Y max ) systematically smaller than the ones 
produced by both versions of SIBYLL. This is due to the 
much higher average particle multiplicity generated by 
QGSjet98 and the lower elasticity compared to SIBYLL. 
These two features are responsible for the accelerated 
shower development in QGSjet98. An interesting feature 
is that the proton- proton cross section in SIBYLL 2.1 is 
~ 25% larger than the one predicted by QGSjet98 at 10 20 
eV, however the larger multiplicity and smaller elasticity 
of the latter still dominate producing a smaller (Y max ) . 
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Finally the mean energy fraction carried by 7r°'s, e s 
and photons is shown in Fig. |8[ Interestingly both models 
predict that the same fraction of the projectile energy 
is transferred to the electromagnetic shower component 
at high energy. However, the electromagnetic showers 
in a SIBYLL 2.1 simulation are more energetic and less 
numerous than with QGSjet98. 



FIG. 9: Average depth of maximum (X max ) of proton showers 
as a function of primary energy. The lines represent 5,000 
events generated by our one dimensional method, at 9 = 45°, 
using SIBYLL 1.7 (dotted), SIBYLL 2.1 (solid) and QGSjet98 
(dashed). The symbols show the values of (X max ) averaged 
over 500 showers obtained with CORSIKA using the thinning 
procedure. 

The width of the Y max distribution is a measure of 
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FIG. 10: Fluctuation of the position of the shower maximum, 
a = yj (^iLix) — (Y max ) 2 . The curves have the same meaning 
as in Fig. 



the fluctuations of the position of the shower maximum. 
As shown in Fig. [h], the fluctuations become less im- 
portant at very high energy. First of all, the fluctua- 
tions due to the position of the first interaction point 
are smaller at high energy due to the large cross section 
(small mean free path). Secondly, the large multiplicity 
of secondary particles produces a correspondingly larger 
number of subshowers. Individual subshowers will show 
considerable profile fluctuations as observed at lower en- 
ergy, however, due to their large number the total shower 
profile exhibits much smaller fluctuations. 

We have verified that the LPM effect Q ||, ||] 
doesn't affect (Y max ) for proton energies below 3 x lfP° 
eV in agreement with |Q . The values of ( Y max ) in pro- 
ton showers at energy 3 x 10 20 and in proton showers at 
the same energy but with the LPM artificially "turned 
off" are equal within ~ 1%. The large multiplicity of 
hadronic interactions at energies above the scale at which 
the LPM is important, is largely responsible for this small 
difference, because it reduces the energy of the neutral 
pions whose decays are the dominant channel for pro- 
duction of high energy photons in the shower. Neutral 
pions then do not produce energetic enough photons as to 
show strong LPM characteristics. Even if a high energy 
neutral pion is created, for instance in diffractive interac- 
tions in which the multiplicity is low, at energies above 
~ 10 19 eV interactions of neutral pions dominate over 
decay and hence the production of high energy photons 
is suppressed Q . 

Numerical values of (X max ) and (5 max ) are presented 
in Table |l| for vertical proton induced showers. A com- 
parison between ( Y max ) in this table and in Fig. ^ reveals 
its weak dependence on the zenith angle in the angu- 
lar range 8 = 0° — 45°. (A max ) is fairly insensitive to 



changes in atmospheric density profile from 8 — 0° to 
8 = 45° and hence it is approximately the same when 
expressed in g/cm 2 . {S max } also shows a weak depen- 
dence on the zenith angle and it is remarkably indepen- 
dent of the hadronic interaction models adopted in this 
work. SIBYLL 2.1 produces (S mayi ) values smaller than 
those predicted by SIBYLL 1.7 by 1% in the whole energy 
range shown in the table. An interesting aspect about 
the behavior of S max /E with primary energy is that it 
increases up to energies of ~ 10 17 eV and decreases after 
that for all three models. 
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FIG. 11: Distribution of S max normalized by the primary 
energy in GeV. Results are shown for 5,000 primary proton 
showers of energies 10 18 eV (bottom panel) and 10 20 eV (top 
panel), with zenith angle 6 = 0° calculated with the hybrid 
method using SIBYLL 1.7 (dotted), SIBYLL 2.1 (solid), and 
QGSjet98 (dashed). 

In Figs. | and we compare our predictions for proton 
showers to those obtained in the framework of the COR- 
SIKA code using similar (or identical) hadronic interac- 
tion models [§2 Each of the points generated 
with CORSIKA in Figs. § and |o|, represents the mean 
value of Y max over 500 showers using the thinning pro- 
cedure. The values of (Y max ) and a calculated by both 
codes for the same models are in very good agreement 
p9| , within the larger statistical uncertainty of this par- 
ticular CORSIKA calculation. This provides us a further 
check on the validity of the hybrid simulation method. 
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TABLE II: Mean depth of shower maximum development, (X max ), and shower size at depth of maximum (S m ax), in proton- 
initiated shower with incident zenith angle 9 — 0°. Each energy represents 5,000 showers simulated with the hybrid method 
using SIBYLL 1.7, SIBYLL 2.1, and QGSjet98. The width of the corresponding distributions is given in parenthesis. 



Model 


SIBYLL 1.7 


SIBYLL2.1 


QGSjet98 


lg(E/eV) 


<A max ) [g/cm 2 ; 


(S ma x)/£ [GeV- 




<Y max )[g/cm 2 ; 


{S m )/E [GeV" 




<X max }[g/cm 2 ] (S max )/£ [GeV" 




14.0 


530 (101) 


0.691 (1.10 x 10" 




507 (96) 


0.688 (1.03 x 10" 




499 (95) 


0.685 (9.93 x 10" 


- 2 ) 


15.0 


592 (86) 


0.719 (8.44 x 10" 


- 2 ) 


571 (82) 


0.719 (7.63 x 10" 


- 2 ) 


565 (86) 


0.724 (7.06 x 10" 


- 2 ) 


16.0 


647 (72) 


0.735 (6.07 x 10" 


- 2 ) 


626 (71) 


0.734 (5.57 x 10" 


- 2 ) 


625 (78) 


0.736 (5.25 x 10" 


- 2 ) 


17.0 


706 (64) 


0.739 (4.33 x 10" 


- 2 ) 


684 (64) 


0.737 (4.04 x 10" 


- 2 ) 


677 (70) 


0.738 (3.75 x 10" 


- 2 ) 


18.0 


760 (57) 


0.737 (3.10 x 10" 




740 (58) 


0.734 (2.97 x 10" 




730 (66) 


0.728 (2.85 x 10" 




19.0 


822 (55) 


0.723 (2.80 x 10" 




799 (55) 


0.718 (2.56 x 10" 


- 2 ) 


785 (66) 


0.708 (2.36 x 10" 


- 2 ) 


20.0 


878 (51) 


0.694 (3.46 x 10" 




856 (51) 


0.690 (3.43 x 10" 


- 2 ) 


832 (62) 


0.683 (2.85 x 10" 


- 2 ) 


20.5 


901 (46) 


0.671 (4.29 x 10" 




880 (47) 


0.666 (3.98 x 10" 




853 (56) 


0.662 (3.55 x 10" 


- 2 ) 
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FIG. 12: Distribution of X max . Results are shown for 5,000 
vertical showers generated by primary protons of energies 10 18 
eV (bottom panel) and 10 20 eV (top panel) calculated with 
the hybrid method using SIBYLL 1.7 (dotted), SIBYLL 2.1 
(solid), and QGSjet98 (dashed). 



Fig. [n] shows the distribution of S max normalized by 
the primary energy in GeV. The top (bottom) histogram 
represents 5,000 proton-induced vertical showers at 1 
EeV (100 EeV), calculated using SIBYLL 1.7, SIBYLL 



2.1 and QGSjet98. The numerical values of (S mlix )/E 
are shown in Table ||. The distribution of S max /E is 
clearly not symmetric around the most likely value. This 
is a common feature of the three models and reflects the 
asymmetric fluctuations of the various interaction points 
and secondary particle multiplicities. 

Fig. |l2| shows the distribution of X max , for the same 
shower initial parameters as in Fig. O The distribution 
has an asymmetric shape with a long tail at large values 
of X max . At both energies the tendency of QGSjet98 to 
predict lower values of X max is clearly visible. The differ- 
ence is more apparent when compared to the distribution 
obtained for SIBYLL 1.7. The distribution also reflects 
graphically the larger fluctuations expected for QGSjet98 
compared to SIBYLL. 

The fluctuations in AT max are directly related to the 
relative fraction of diffraction dissociation events gener- 
ated in SIBYLL and QGSjet. In particular showers which 
develop very deep in the atmosphere are typically those 
with a diffractive first interaction. Inelastic diffraction in 
proton-air collisions can be subdivided into coherent and 
incoherent diffraction. The latter process corresponds to 
the interaction of the projectile with a single nucleon of 
the target nucleus and is therefore completely analogous 
to diffraction in proton-proton collisions. As a multi- 
channel eikonal model (52) SIBYLL 2.1 predicts a growth 
of the cross section for diffraction dissociation in proton- 
proton collisions like In s which means that the fraction of 
low- multiplicity events decreases at high energy as 1/ In s 
f33|] . In contrast, in QGSjet the fraction of diffractive 
events is essentially energy independent (more precisely 
proportional to the ratio of the elastic and inelastic cross 
sections) because it is based on the quasi-eikonal approx- 
imation |54| . From theoretical grounds the quasi-eikonal 
approximation is expected to overestimate the diffractive 
cross section at high energy as it does not implement the 
black disk limit (for a discussion of the black disk limit 
see, for example, |)3j). On the other hand QGSjet ac- 



11 



counts also for coherent diffraction which is neglected in 
SIBYLL. 



C. Elongation Rate 



The elongation rate is defined as |22], 23 

d(X max ) 



D 



10 



dlgE 



(11) 



It describes the energy-dependence of the position of the 
shower maximum. The elongation rate reflects changes in 
the cosmic ray composition as well as features of hadronic 
interaction at high energy. Our interest here is in the re- 
lation between elongation rate and hadronic interactions. 

Most of the charged particles in the shower are elec- 
trons and positrons with energies near the critical en- 
ergy (81 MeV in air) from electromagnetic subshowers 
initiated by photons from 7r°-decay. The mean depth of 
maximum for an electromagnetic shower initiated by a 

35fl 



photon with energy Ey is | 



A ln£L + C 



(12) 



where Xq ps 37g/cm 2 is the radiation length in air. The 
elongation rate for an electromagnetic shower is thus 
Dfg 1 = ln(10) x X Q « 85 g/cm 2 . 

A proton-initiated shower consists of a hadronic 
core feeding the electromagnetic component primarily 
through 7T° production. In the approximation of a 
hadronic interaction model that obeys Feynman scaling 
with energy-independent cross sections, the energy split- 
ting in the hadronic skeleton of the shower is independent 
of energy (i.e. it scales with energy). As a consequence, 
since the electromagnetic component is dominated by the 
earliest (i.e. most energetic) generations of hadronic in- 
teractions, under these assumptions the elongation rate 
of the hadronic shower is also -Df™. In general, for an in- 
cident nucleus of mass A and total energy Eq (including 
protons with A = 1) the depth of maximum is 



(X n 



X ln(E /A) + X A , 



(13) 



where is the interaction length of the primary particle. 
If the composition changes with energy, then (A) depends 
on energy and the elongation rate changes accordingly. 

In qualitative analyses of the role of hadronic interac- 
tions in air shower development, an approach analogous 
to the treatment of nuclei has often been used. The depth 
of maximum for a proton shower is expressed as 



(X^(E)) = (X^JE/(n)))+X N , 



(14) 



where (n) is related to the multiplicity of secondaries in 
the high-energy hadronic interactions in the cascade. The 
situation is, however, essentially more complicated than 
for a primary nucleus in which the energy is to a good 
approximation simply divided into A equal parts. In a 
hadronic cascade instead there is a hierarchy of energies 



of secondary particles in each interaction, and a similar 
(approximately geometric) hierarchy of interaction ener- 
gies in the cascade. In this case (n) has to be understood 
as some kind of "effective" multiplicity, which does not 
have a straightforward definition in general. 
The elongation rate derived from Eq. ([l4|) is 



d(X^(E)) 
d\gE 



ln(10)A 



dln( 



d\nE 



dX 



N 



dlgE 1 



(15) 



which corresponds to the form given by Linsley and Wat- 
son M, 



D W =H10)X (1-B n -B x ), 



with 



B n = 



(iln(n) 
dlnE 



B x = - 



Ajy din Ajv 
Aq d\nE ' 



(16) 



(17) 



For a hadronic interaction model with a multiplicity 
dependence of (n) = uqE s one gets B n = S provided 
all secondaries having the same energy, which is not the 
case. 
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FIG. 13: Elongation rates, d(X ma x)/cilg E, calculated numer- 
ically using showers simulated with the hybrid method (see 
text). 



Because Eq. ( |16| ) is often used to estimate the elon- 
gation rate (see, for example, [p6[), it is worthwhile to 
compare our results with this parametrization. Fig. [l3| 
shows the elongation rate for SIBYLL and QGSjet show- 
ers as derived from the detailed shower simulation. All 
the models show an initial decline from the low-energy 
scaling regime as expected. Then, above 10 15 eV the 
elongation rate for SIBYLL is nearly constant while that 
for QGSjet continues to decline. In addition SIBYLL 
1.7 has a sharp drop of the elongation rate at ultra-high 
energy, which we explain below. In contrast, if we dif- 
ferentiate the curves in Figs. [| and [7] for cross section 
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and multiplicity and calculate the elongation rate from 
(|l6|), using the assumption of equal sharing of energy 
among the secondaries, we get completely misleading re- 
sults, particularly at low energy. For example, for QGSjet 
the predicted elongation rate is about 60 g/cm 2 over the 
entire energy range. The situation is similar for showers 
simulated with SIBYLL 2.1. 

The use of the total particle multiplicity for (n) is not 
so bad at high energy because the scaling violation is 
fully developed in the high energy part of the shower. 
However, it is important to note that in Eq. (16) the vi- 
olation of Fcynman scaling and the energy-dependence 
of the cross sections are taken into account only for the 
first interaction. All subsequent hadronic interactions 
are assumed to be characterized by Feynman scaling and 
constant interaction cross sections. Thus Eq. (16) is ex- 
pected to be a good approximation only in an interme- 
diate energy range around 10 15 — 10 16 eV. 

At higher energy the energy-dependence of the subse- 
quent hadronic interactions becomes important. As an 
illustration, we consider a toy model in which all final 
state particles of the first proton interaction in air are 
charged pions and have the same energy. At high en- 
ergy all pions will interact before decaying. As a first 
approximation we can write 



(X^(E)) = (X^(E/(n))) + X N 



where now the position of maximum of pion induced sec- 
ondary showers is written at r.h.s. Using Eq. ( [l5| ) to de- 
scribe the pion showers one gets for the elongation rate 
of the entire shower 



d\gE 
ln(10)X 



d\n(n{E)) d\n(n(E/n(E))) 



dlnE 



d\uE 



+ 



dX„ 



d\gE d\gE 



In a model with a power-law increase of the multiplicity 
with index 8 this simplifies to 



d\gE 



ln(10)X [1 - 25 + 5 2 ] 



dX 



N 



d\„ 



d\gE d\gE 
(20) 



Usingagain the cross sections and multiplicities shown in 
Figs. g| and (?] one gets elongation rates of about 43 and 56 
g/cm for QGSjet and SIBYLL, respectively. Given the 
simplicity of the model the predictions are remarkably 
close to the results of the full simulation above 10 19 eV. 
This could be the result of a cancellation of two effects: 
on one hand only two successive hadronic interactions 
were assumed to be energy-dependent and, on the other 
hand, the scaling violation in these interactions was over- 
estimated by using the total particle multiplicity in jl9| ) 
and the uniform energy sharing. 



Finally it should be mentioned that the sudden drop 
of the elongation rate of the showers simulated with 
SIBYLL 1.7 is due to the onset of the interaction of 
neutral pions. At energies above 10 195 eV a substan- 
tial number of 7r°'s does not decay but interacts because 
of the enormous Lorentz dilation. This effect reduces the 
mean energy of the particles which feed the electromag- 
netic component of the shower. The change in elongation 
rate is most prominent in SIBYLL 1.7 because it gener- 
ates more fast (interacting) neutral pions. 



D. Number of Muons 

The number of muons in a shower is an important ob- 
servable which depends strongly on the mass of the pri- 
mary particle and is used in the studies of the elemental 
composition of cosmic rays. It also directly reflects the 
hadronic component of the shower and hence it is a sen- 
sitive probe of the hadronic interactions. 

We have calculated the average number of muons at 
sea level ((A^)) with energies above £* hr =0.3, 1, 3, 10 
and 30 GeV, in proton-initiated showers at zenith angle 
9 = 0° (0 = 45°) for the hadronic models SIBYLL 1.7, 
SIBYLL 2.1, and QGSjet98. Fig. || shows the energy 
dependence of the average number of muons normalized 



(18) to the primary energy for the three models. (A/,) follows 
approximately a simple power law (E/E c ) a for energies 
above ~ 10 14 eV. 

This can be understood on the basis of Hcitler's model 



]57[ (see also the discussion in 58 ) by assuming that 
each hadronic interaction produces in average (n to t) sec- 
ondaries of approximately the same energy. The multi- 
plication of the number of charged pions in a shower con- 
tinues until the pions reach a critical energy, E c , at which 
they are assumed to decay. After N generations (i.e. sub- 
sequent interactions) the energy of the pions reaches the 
critical energy E c = E/{n% Q t) N . The number of muons 



(19) from decaying charged pions is thus A M 



inating N gives 



N„ 



ln(n„ 



ln(ntot) ' 



Elim- 



(21) 



which is the well-known power-law found in data. The in- 
dex a can be calculated by using (n 7T ±) ss §("tot)j which 
gives values for a in the range from 0.85 to 0.92. (Assum- 
ing that the charged pion multiplicity is less than 2/3 of 
the total multiplicity decreases the values predicted for 

Over the entire energy range from 10 14 eV to more 
than 10 20 eV a single power law parametrization can be 
used to describe the muon multiplicities for all energy 
threshold considered here. In Tabs. Ill and IV we show 



the corresponding fit parameters for showers of and 45° 
zenith angle, respectively. As expected the critical energy 
increases with the muon threshold energy. The energy- 
dependence of the muon multiplicity is the steepest for 
low-energy muons. For a given muon energy threshold, 
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TABLE III: Parameters a and E c obtained by fitting the number of muons in vertical showers at sea level using a power law of 
the form jV M = (E/E c ) a . The numerical values of the parameters are presented for the three hadronic models and for muons 
with energy above 0.3, 1, 3, 10 and 30 GeV. 



Model 


SIBYLL 1.7 


SIBYLL 2.1 


QGSjet98 


E^ [GeV] 


0.3 1 3 10 30 


0.3 1 3 10 30 


0.3 1 3 10 30 


a 


0.886 0.877 0.869 0.857 0.846 


0.901 0.893 0.884 0.872 0.861 


0.920 0.913 0.904 0.893 0.882 


E c [GeV] 


35 43 67 162 594 


39 47 70 161 555 


44 53 79 182 638 



TABLE IV: Parameters a and E c obtained by fitting the number of muons in showers with 6 = 45° at sea level using a power 
law of the form JV M = (E/E c ) a . The numerical values of the parameters are presented for the three hadronic models and for 
muons with energy above 0.3, 1, 3, 10 and 30 GeV. 



Model 


SIBYLL 1.7 


SIBYLL 2.1 


QGSjet98 


Ef [GeV] 


0.3 1 3 10 30 


0.3 1 3 10 30 


0.3 1 3 10 30 


a 


0.891 0.886 0.877 0.867 0.853 


0.902 0.897 0.890 0.878 0.865 


0.921 0.916 0.909 0.899 0.887 


E c [GeV] 


65 72 94 195 562 


67 74 96 185 523 


77 83 107 209 607 



the numerical value of a tends to be the highest for the 
QGSjet98 model. 

Already from the simple model discussed above it 
is clear that the power-law index should be energy- 
dependent because the multiplicity of the secondary par- 
ticles increases with energy. Indeed, a careful inspection 
of the energy dependence of (iV M ) shows that the power 
law index a increases with the primary energy. However, 
the observed relative deviation from a single power law 
is always less than 15%. It is the regularity of this devi- 
ation and the aforementioned physics motivation which 
makes it worthwhile to consider the following, alternative 
parametrization. 

The power-law index is taken to be energy-dependent 
with 



a(E) 



ln(3/2) 
ln(ra c ff) 



(22) 



where n e g is the geometric average of the charged pion 
multiplicity of N successive hadronic interactions. By 
construction this effective multiplicity has a weak energy 
dependence, which we approximate by 



ln(n e ff) ~ no + n\ In 



E 



E a = 10 14 eV. (23) 



To make the numerical values of a(E) more transparent 
we express the parameters no and n\ in terms of power- 
law indices ao = ol{Eq) and a,\ = a(E\ = 10 20 eV) 



n = 



Of) 



1 - a 

ln(3/2) 
HEhjEo) 



ln(3/2) 



OLQ 



1 — OL\ 1 — OL0 



(24) 
(25) 



This alternative muon multiplicity parametrization has 
only three free parameters, the indices ao^ ol\ and the 



critical energy E c . It gives considerably better fits to the 
simulation data than the single power-law parametriza- 
tion ((|^). The numerical values obtained by fitting the 
output of the hybrid simulations are shown in Tables |v| 
and VI. The relative uncertainties of the parameters ao. 



a\ are about 1% and 10 to 15% for E c . 

The QGSjet98 model shows the biggest change of the 
power law index from ao to a±. Muon production in 
SIBYLL 2.1 is the closest to a simple power law. The 
general trend for all three models is that the power law 
index decreases with the muon threshold energy. 

The absolute number of muons differs from model to 
model. SIBYLL 2.1 produces more muons than SIBYLL 
1.7 but still less than QGSjet98 at all energies. The dif- 
ferences between the three models increase with energy 
and reach maximum at 10 20 eV. Table VII gives the ra- 
tios of generated by SIBYLL 1.7 and QGSjet98 at 
sea level to those generated by SIBYLL 2.1 in vertical 
showers at primary energies of 10 15 eV and 10 20 eV. 

It is interesting to observe the dependence of these 
differences on the muon threshold energy. While for 
SIBYLL 1.7 the ratio decreases monotonically with the 
threshold energy, the QGSjet98/SIBYLL 2.1 ratio shows 
a more complex behavior. The enhanced production of 
low energy muons in QGSjet98 is related to the higher 
charged multiplicity of the model in the 100 - 1000 GeV 
range. The differences between the two models decrease 
for E^ r of 30 GeV. 

The number of muons at sea level is sensitive to the 
incident zenith angle. Two competing processes - muon 
production and muon energy loss and decay - determine 
the dependence on zenith angle. With increasing zenith 
angle both the grammage in which showers develop and 
the distance to the observation level increase. Some ad- 
ditional muons are generated in inclined showers due to 
the larger number of interactions, but also a large fraction 
of the low energy muons (below ~ 3 GeV) decay before 
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TABLE V: Parameters «o, o?i and E c obtained by fitting the number of muons in vertical showers at sea level to Eq. (|22j). The 
numerical values of the parameters are presented for the three hadronic models and for muons with energy above 0.3, 1, 3, 10 
and 30 GeV. 



Model 


SIBYLL 1.7 


SIBYLL 2.1 


QGSjet98 


Ef 1 [GeV] 


0.3 1 3 10 30 


0.3 1 3 10 30 


0.3 1 3 10 30 


ato 


0.858 0.838 0.819 0.780 0.745 


0.887 0.870 0.850 0.820 0.787 


0.855 0.834 0.809 0.775 0.736 


ai 


0.874 0.861 0.849 0.827 0.809 


0.895 0.883 0.870 0.852 0.834 


0.892 0.879 0.864 0.846 0.828 


E c [GeV] 


26 28 39 74 238 


33 36 49 97 291 


22 23 29 57 179 



TABLE VI: Parameters eto, ai and E c obtained by fitting the number muons in inclined showers (9 — 45°) at sea level to 
Eq. (p2j). The numerical values of the parameters are presented for the three hadronic models and for muons with energy above 
0.3, 1, 3, 10 and 30 GeV. 



Model 


SIBYLL 1.7 


SIBYLL 2.1 


QGSjet98 


Eft" [GeV] 


0.3 1 3 10 30 


0.3 1 3 10 30 


0.3 1 3 10 30 


cto 


0.873 0.863 0.842 0.819 0.764 


0.892 0.883 0.867 0.839 0.802 


0.858 0.848 0.828 0.794 0.751 


Ql 


0.884 0.876 0.863 0.849 0.821 


0.898 0.891 0.881 0.863 0.842 


0.895 0.888 0.876 0.858 0.837 


E c [GeV] 


54 57 66 124 255 


61 64 76 128 304 


41 42 48 78 189 



TABLE VII: Ratios of (JV^,) at sea level generated in vertical 
showers by SIBYLL 1.7 and QGSjet98 to those generated by 
SIBYLL 2.1 (= 1). 



E [eV] 


10 15 


10 20 


Ef 1 [GeV] 


0.3 3 30 


0.3 3 30 


SIBYLL 1.7 
QGSjet98 


0.96 0.92 0.87 
1.11 1.12 1.06 


0.83 0.81 0.78 
1.37 1.41 1.35 



reaching sea level. Decays win the competition and the 
number of low energy muons decreases with zenith angle. 
At energies above ~ 10 GeV, however, most of the muons 
cross the whole atmosphere without decaying, and their 
number at sea level is less sensitive to the injection angle. 

This behavior is illustrated in Fig. [l5] which shows the 
distribution of the number of muons at sea level for differ- 
ent E^ and zenith angles of 0° and 45° . Each histogram 
represents 5,000 showers initiated by primary protons at 
1 EeV using the SIBYLL 2.1 model. 

At energy above 30 GeV practically all muons cross 
the atmosphere without decaying. The difference in the 
number of muons above 30 GeV between the two zenith 
angles, depicted in Fig. is then determined by muon 
production. At large zenith angles shower particles travel 
for a longer time in a more tenuous atmosphere and hence 
the charged pions have a smaller probability of interac- 
tion. As a result more muons are produced at = 45° 
than at 6 = 0°. 

SIBYLL 2.1 and QGSjet98 predict similar fluctuations 
in the number of muons. At E = 10 18 eV the width of the 
shower distribution in muons obtained with QGSjet98 is 



only ~ 7% larger than in SIBYLL 2.1 for all muon energy 
thresholds. The difference in the widths at 10 18 eV as 
obtained with QGSjet98 and SIBYLL 1.7 is larger and 
increases from - 17% at Ef 1 = 0.3 GeV to - 27% at 
Ef 1 = 30 GeV. 

IV. SUMMARY AND OUTLOOK 

We have presented an efficient, one-dimensional hy- 
brid method to simulate the development of extensive 
air showers. The combination of Monte Carlo techniques 
for the interactions of the shower particles above a cer- 
tain hybrid energy threshold with a presimulated library 
of pion-induced showers, allows us to simulate the devel- 
opment of large statistical samples of air showers up to 
the highest energies observed. 

Previously developed hybrid methods use the average 
longitudinal development to describe the numerous sub- 
threshold showers and are usually limited to the calcula- 
tion of the total number of electromagnetic particles. In 
this paper we have presented a method that accounts 
for fluctuations in the shower development as well as 
the correlations between the different parameters describ- 
ing the electromagnetic and muon components of EAS. 
By comparing direct simulations with hybrid-simulated 
showers we have determined that the correlation between 
the hadronic and the electromagnetic component is also 
well reproduced with our method. In particular the hy- 
brid method correctly describes the correlation between 
the number of muons and the shower size at observation 
level, which is of special relevance to studies of the cosmic 
ray composition. 
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FIG. 14: Average number of muons at sea level (iV M ), obtained 
in proton showers with zenith angle 9 = 0°. Each energy 
represents 5,000 showers simulated with the hybrid method. 
The solid (dotted) line represents the values obtained with 
SIBYLL 2.1 (SIBYLL 1.7), while the dashed line illustrates 
the values for QGSjet98. Panels (a), (b) and (c) show the 
average number of muons with energy above 30 GeV, 3 GeV 
and 0.3 GeV respectively. 



We have studied the influence of different hadronic in- 
teraction models, namely SIBYLL 1.7, SIBYLL 2.1 and 
QGSjet98, on shower observables which are relevant for 
the determination of the energy and chemical composi- 
tion of the primary cosmic ray flux. We have presented 
average values of A" max , S m & x and the number of muons 
above 0.3, 1, 3, 10 and 30 GeV at sea level, as well as the 
fluctuations of these quantities. The mean muon multi- 
plicities were analyzed with two different models: 

(i) a simple power-law parametrization, which describes 
the simulation results with a relative accuracy of better 
than 10% (15% for £* hr = 30 GeV), and 

(ii) a model with a slowly changing power-law index, 
which gives an excellent description of the data. 

The relation between the features of the interaction 
models and the shower observables has been extensively 
discussed. We stress the influence of the different extrap- 
olations of the hadronic models to the highest energies on 
the features of the electromagnetic and hadronic compo- 
nent of the shower, and the influence of the differences 
between the models on the number of muons predicted 
by them. Some of these differences exist already at low 
energies and affect the average numbers of low energy 
muons. 

In QCD-inspired models such as SIBYLL and QGSjet 
the predictions on cross sections are inherently linked to 
the size of Feynman scaling violation, and hence multi- 




le+06 



4e+06 
N„ 
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FIG. 15: Shower distribution in number of muons at sea level. 
The results are obtained for 5,000 primary proton showers of 
energy 10 18 eV for different muon energy thresholds. The 
solid line represents vertical showers, while the dotted line 
illustrates showers with zenith angle — 45°. All showers 
were simulated using SIBYLL 2.1. 



plicity, implemented in the model. A model with a steep 
energy-dependence of the hadron-air cross section is usu- 
ally characterized by a moderate increase of the multi- 
plicity. Concerning the position of the shower maximum 
the effect of large Feynman scaling violation (or a steeply 
rising multiplicity) is similar to that of a steeply rising 
cross section. This is the reason why the (X max ) predic- 
tions of SIBYLL 2.1 and QGSjet98 are rather similar over 
a wide energy range. On the other hand, the number of 
muons at sea level reflects the multiplicity of low-energy 
hadrons produced in a shower but depends only weakly 
on the hadronic cross sections. Therefore, showers simu- 
lated with QGSjet produce consistently more low-energy 
muons than SIBYLL showers. 

As another application of our method wc have studied 
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the influence of the multiplicity, inelasticity, proton-air 
cross section on the elongation rate of proton-initiated 
showers. We find that the elongation rate has a complex 
dependence on the scaling violation and cross section be- 
haviour of hadronic interaction models. Again, a steeply 
rising cross section leads to a decrease of the elongation 
rate qualitatively similar to a steeply rising multiplicity. 
Furthermore, a threshold-like behaviour is observed at 
extremely high energy. The onset of the hadronic inter- 
action of neutral pions, which always decay at low energy, 
leads to a significant decrease of the elongation rate. 

In forthcoming work we will apply our hybrid method 
to the determination of the proton-air cross section in 
experiments that are able to measure the muon and elec- 
tromagnetic components at fixed depth, as well as in ex- 
periments capable of measuring the distribution of X max . 
Furthermore, we will exploit the fastness of our method 
to simulate large statistical samples of showers initiated 
by heavy nuclei, with the aim to predict observables that 



help studying the composition of the cosmic ray flux. 
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